iSSF Vegetation Model Plotting

Here we are plotting the results of integrated step selection function (iSSF) models fitted to water buffalo in Arnhem Land. These models contain vegetation categories as a covariate, with interactions for movement parameters.

Load packages

library(tidyverse)
library(TwoStepCLogit) # to make population estimates
library(sjPlot) # to plot model estimates
library(glmmTMB)# for the mixed models 
library(tictoc)
library(beepr)
library(lemon)
library(circular)
library(amt)

Import data

ssf_dat <- read_csv("pheno_start_end_dat_ssf10rs_2024-03-14.csv")

New names:
Rows: 1062908 Columns: 42
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," chr
(4): pheno_start, pheno_end, season, day_period dbl (35): ...1, id, burst_,
x1_, x2_, y1_, y2_, sl_, ta_, dt_, step_id_, sr... lgl (1): case_ dttm (2):
t1_, t2_
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`

# ssf_dat <- read_csv("pheno_forest_start_end_dat_ssf10rs_2025-03-02.csv")

Prepare the data

# scale the covariates
ssf_dat <- ssf_dat %>% mutate(srtm_start_scaled = scale(srtm_start),
                            srtm_end_scaled = scale(srtm_end),
                            water_dist_start_scaled = scale(water_dist_start),
                            water_dist_end_scaled = scale(water_dist_end),
                            ndvi_start_scaled = scale(ndvi_start),
                            ndvi_end_scaled = scale(ndvi_end),
                            nbr_start_scaled = scale(nbr_start),
                            nbr_end_scaled = scale(nbr_end),
                            sl_scaled = scale(sl_),
                            log_sl_scaled = scale(log_sl_),
                            cos_ta_scaled = scale(cos_ta_)
)

# get scaling factors
sl_scale_sd <- attr(ssf_dat$sl_scaled, "scaled:scale")
log_sl_scale_sd <- attr(ssf_dat$log_sl_scaled, "scaled:scale")
cos_ta_scale_sd <- attr(ssf_dat$cos_ta_scaled, "scaled:scale")

# rename the pheno categories to set a different reference level
ssf_dat$pheno_end <- relevel(as.factor(ssf_dat$pheno_end), ref = "shrub_grass")
ssf_dat$pheno_start <- relevel(as.factor(ssf_dat$pheno_start), ref = "shrub_grass")

unique(ssf_dat$pheno_start)

[1] shrub_grass   dry_grassland tree_grass    bare_ground   open_forest  
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grass

unique(ssf_dat$pheno_end)

[1] shrub_grass   dry_grassland tree_grass    bare_ground   open_forest  
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grass

ssf_dat <- ssf_dat %>% filter(!pheno_start == "saltpan_mudlfat")

# separate the buffalo by season 
ssfdat_wet <- ssf_dat %>% filter(season=="wet") # wet season
ssfdat_dry <- ssf_dat %>% filter(season=="dry") # dry season only 

Plot the used-available distributions

Dry season

# checking for all individuals
buffalo_ids <- unique(ssf_dat$id)

# inspect the data for dry season 
ssfdat_dry %>% 
  group_by(case_, pheno_end) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Dry season", subtitle = "All buffalo")+
  scale_fill_brewer(palette = "Paired", name="Case", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

ggsave("outputs/plots/used-available_all_dry.png", device = "png", 
       width = 150, height = 90, units = "mm", scale = 1.5, dpi = 1000)

Inspect the data for each buffalo

for(i in 1:length(buffalo_ids)) {
  print(ssfdat_dry %>% filter(id == buffalo_ids[i]) %>%
          group_by(case_, pheno_end) %>% 
          summarize(n = n()) %>% 
          mutate(prop = n / sum(n), 
                 label = paste0(round(prop * 100, 1), "%")) %>% 
          ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
          geom_col(position = position_dodge2()) +
          geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
          labs(x = "Land use class", y = "Proportion", fill = "case_")+
          ggtitle("Dry season", subtitle = paste0("Buffalo ", buffalo_ids[i])) +
          scale_fill_brewer(palette = "Paired", name="Case", 
                            breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
          theme_light() +
          theme(axis.text.x = element_text(angle = 20, hjust = 1))
  )
}

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Wet season

ssfdat_wet %>% 
  group_by(case_, pheno_end) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "case_")+
  ggtitle("Wet season", subtitle = "All buffalo")+
  scale_fill_brewer(palette = "Paired", name="Case", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

ggsave("outputs/plots/used-available_all_wet.png", device = "png", 
       width = 150, height = 90, units = "mm", scale = 1.5, dpi = 1000)

For each individual

for(i in 1:length(buffalo_ids)) {
  print(ssfdat_wet %>% filter(id == buffalo_ids[i]) %>%
          group_by(case_, pheno_end) %>% 
          summarize(n = n()) %>% 
          mutate(prop = n / sum(n), 
                 label = paste0(round(prop * 100, 1), "%")) %>% 
          ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
          geom_col(position = position_dodge2()) +
          geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
          labs(x = "Land use class", y = "Proportion", fill = "case_")+
          ggtitle("Wet season", subtitle = paste0("Buffalo ", buffalo_ids[i]))+
          scale_fill_brewer(palette = "Paired", name="Case", 
                            breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
          theme_light() +
          theme(axis.text.x = element_text(angle = 20, hjust = 1))
  )
}

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Read in the saved model objects

buffalo.pheno.dry.mvmt <- readRDS("buffalo_pheno_mixed_corr_mvmt_10rs_shrubgrassref_dry_2025-05-27.rds")
buffalo.pheno.wet.mvmt <- readRDS("buffalo_pheno_mixed_corr_mvmt_10rs_shrubgrassref_wet_2025-05-27.rds")

Model summaries

# Summaries are also printed below
# summary(buffalo.pheno.dry.mvmt)
# summary(buffalo.pheno.wet.mvmt)

Model fitting - interacting with movement

Dry season

tic()
buffalo.tmp <- glmmTMB(case_ ~ 
                         pheno_end + 
                         sl_scaled +
                         log_sl_scaled +
                         cos_ta_scaled +
                         pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
                         (1 | step_id) + 
                         (0 + pheno_end | id),
                         # diag(0 + pheno_end | id),
                       family = poisson(), 
                       data = ssfdat_dry, 
                       doFit=FALSE,
                       control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))

#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA,1:(nvarparm-1))))

# to render the script we comment out the model fitting, as we are loading a saved model
# buffalo.pheno.dry.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)

summary(buffalo.pheno.dry.mvmt) 

 Family: poisson  ( log )
Formula:          
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +  
    (1 | step_id) + (0 + pheno_end | id)
Data: ssfdat_dry

      AIC       BIC    logLik -2*log(L)  df.resid 
       NA        NA        NA        NA    620288 

Random effects:

Conditional model:
 Groups  Name                   Variance  Std.Dev.  Corr                    
 step_id (Intercept)            1.000e+06 1.000e+03                         
 id      pheno_endshrub_grass   6.061e-02 2.462e-01                         
         pheno_endbare_ground   2.876e-02 1.696e-01  0.46                   
         pheno_enddry_grassland 4.862e-02 2.205e-01  0.47  0.43             
         pheno_endopen_forest   1.480e-01 3.847e-01  0.03 -0.13 -0.26       
         pheno_endtree_grass    4.412e-03 6.642e-02 -0.68 -0.18  0.17 -0.43 
Number of obs: 620323, groups:  step_id, 56393; id, 15

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)
(Intercept)                            -3.037772   4.212555  -0.721  0.47083
pheno_endbare_ground                    0.292761   0.062477   4.686 2.79e-06
pheno_enddry_grassland                  0.030961   0.065966   0.469  0.63882
pheno_endopen_forest                   -0.020758   0.119447  -0.174  0.86204
pheno_endtree_grass                    -0.007088   0.079622  -0.089  0.92907
sl_scaled                               0.017326   0.007989   2.169  0.03011
log_sl_scaled                          -0.003433   0.008243  -0.417  0.67704
cos_ta_scaled                          -0.001176   0.006586  -0.179  0.85828
sl_scaled:pheno_startbare_ground       -0.024474   0.017733  -1.380  0.16754
sl_scaled:pheno_startdry_grassland      0.057859   0.016770   3.450  0.00056
sl_scaled:pheno_startopen_forest       -0.033725   0.016938  -1.991  0.04648
sl_scaled:pheno_starttree_grass         0.018940   0.018619   1.017  0.30905
log_sl_scaled:pheno_startbare_ground   -0.040437   0.015810  -2.558  0.01054
log_sl_scaled:pheno_startdry_grassland -0.093955   0.017168  -5.473 4.43e-08
log_sl_scaled:pheno_startopen_forest    0.074908   0.017827   4.202 2.65e-05
log_sl_scaled:pheno_starttree_grass     0.013399   0.019710   0.680  0.49662
cos_ta_scaled:pheno_startbare_ground   -0.046387   0.012831  -3.615  0.00030
cos_ta_scaled:pheno_startdry_grassland  0.009336   0.013844   0.674  0.50006
cos_ta_scaled:pheno_startopen_forest   -0.037983   0.013482  -2.817  0.00484
cos_ta_scaled:pheno_starttree_grass    -0.035660   0.014887  -2.395  0.01660
                                          
(Intercept)                               
pheno_endbare_ground                   ***
pheno_enddry_grassland                    
pheno_endopen_forest                      
pheno_endtree_grass                       
sl_scaled                              *  
log_sl_scaled                             
cos_ta_scaled                             
sl_scaled:pheno_startbare_ground          
sl_scaled:pheno_startdry_grassland     ***
sl_scaled:pheno_startopen_forest       *  
sl_scaled:pheno_starttree_grass           
log_sl_scaled:pheno_startbare_ground   *  
log_sl_scaled:pheno_startdry_grassland ***
log_sl_scaled:pheno_startopen_forest   ***
log_sl_scaled:pheno_starttree_grass       
cos_ta_scaled:pheno_startbare_ground   ***
cos_ta_scaled:pheno_startdry_grassland    
cos_ta_scaled:pheno_startopen_forest   ** 
cos_ta_scaled:pheno_starttree_grass    *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

toc()

4.288 sec elapsed

beep(sound = 2)

# plot result as incidence rate ratios
# plot_model(buffalo.pheno.dry.mvmt)
plot_model(buffalo.pheno.dry.mvmt, transform = NULL)

# saveRDS(buffalo.pheno.dry.mvmt, file = paste0("buffalo_pheno_mixed_corr_mvmt_10rs_shrubgrassref_dry_", Sys.Date(), ".rds"))
# paste0("Model location: buffalo_pheno_mixed_corr_mvmt_10rs_shrubgrassref_dry_", Sys.Date(), ".rds")

Wet season

tic()
buffalo.tmp <- glmmTMB(case_ ~ 
                         pheno_end + 
                         sl_scaled +
                         log_sl_scaled +
                         cos_ta_scaled +
                         pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
                         (1 | step_id) + 
                         (0 + pheno_end | id),
                         # diag(0 + pheno_end | id),
                       family=poisson(), 
                       data = ssfdat_wet, 
                       doFit=FALSE,
                       control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))

#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA,1:(nvarparm-1))))

# to render the script we comment out the model fitting, as we are loading a saved model
# buffalo.pheno.wet.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)

summary(buffalo.pheno.wet.mvmt)

 Family: poisson  ( log )
Formula:          
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +  
    (1 | step_id) + (0 + pheno_end | id)
Data: ssfdat_wet

      AIC       BIC    logLik -2*log(L)  df.resid 
       NA        NA        NA        NA    442550 

Random effects:

Conditional model:
 Groups  Name                   Variance  Std.Dev.  Corr                    
 step_id (Intercept)            1.000e+06 1.000e+03                         
 id      pheno_endshrub_grass   1.920e-02 1.386e-01                         
         pheno_endbare_ground   8.897e-02 2.983e-01  0.10                   
         pheno_enddry_grassland 7.525e-04 2.743e-02 -0.29  0.79             
         pheno_endopen_forest   6.469e-02 2.543e-01  0.42 -0.45 -0.22       
         pheno_endtree_grass    1.771e-02 1.331e-01 -0.46  0.39  0.53 -0.40 
Number of obs: 442585, groups:  step_id, 40235; id, 15

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)
(Intercept)                            -2.884297   4.985526  -0.579 0.562904
pheno_endbare_ground                    0.513949   0.086359   5.951 2.66e-09
pheno_enddry_grassland                  0.075502   0.045540   1.658 0.097333
pheno_endopen_forest                    0.037996   0.069128   0.550 0.582561
pheno_endtree_grass                     0.058427   0.067305   0.868 0.385342
sl_scaled                              -0.001822   0.008832  -0.206 0.836523
log_sl_scaled                           0.153830   0.011178  13.762  < 2e-16
cos_ta_scaled                           0.084196   0.007907  10.648  < 2e-16
sl_scaled:pheno_startbare_ground       -0.009463   0.018312  -0.517 0.605325
sl_scaled:pheno_startdry_grassland      0.007007   0.021415   0.327 0.743514
sl_scaled:pheno_startopen_forest       -0.090825   0.022800  -3.984 6.79e-05
sl_scaled:pheno_starttree_grass        -0.039944   0.025029  -1.596 0.110508
log_sl_scaled:pheno_startbare_ground   -0.074127   0.019116  -3.878 0.000105
log_sl_scaled:pheno_startdry_grassland -0.050749   0.025079  -2.024 0.043016
log_sl_scaled:pheno_startopen_forest   -0.113883   0.023047  -4.941 7.76e-07
log_sl_scaled:pheno_starttree_grass    -0.156723   0.024816  -6.315 2.70e-10
cos_ta_scaled:pheno_startbare_ground   -0.104247   0.014004  -7.444 9.77e-14
cos_ta_scaled:pheno_startdry_grassland -0.013032   0.017470  -0.746 0.455713
cos_ta_scaled:pheno_startopen_forest   -0.156835   0.017240  -9.097  < 2e-16
cos_ta_scaled:pheno_starttree_grass    -0.141395   0.018060  -7.829 4.90e-15
                                          
(Intercept)                               
pheno_endbare_ground                   ***
pheno_enddry_grassland                 .  
pheno_endopen_forest                      
pheno_endtree_grass                       
sl_scaled                                 
log_sl_scaled                          ***
cos_ta_scaled                          ***
sl_scaled:pheno_startbare_ground          
sl_scaled:pheno_startdry_grassland        
sl_scaled:pheno_startopen_forest       ***
sl_scaled:pheno_starttree_grass           
log_sl_scaled:pheno_startbare_ground   ***
log_sl_scaled:pheno_startdry_grassland *  
log_sl_scaled:pheno_startopen_forest   ***
log_sl_scaled:pheno_starttree_grass    ***
cos_ta_scaled:pheno_startbare_ground   ***
cos_ta_scaled:pheno_startdry_grassland    
cos_ta_scaled:pheno_startopen_forest   ***
cos_ta_scaled:pheno_starttree_grass    ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

toc()

2.914 sec elapsed

beep(sound = 2)

# plot result as incidence rate ratios
# plot_model(buffalo.pheno.wet.mvmt)
plot_model(buffalo.pheno.wet.mvmt, transform = NULL)

# saveRDS(buffalo.pheno.wet.mvmt, file = paste0("buffalo_pheno_mixed_corr_mvmt_10rs_shrubgrassref_wet_", Sys.Date(), ".rds"))
# paste0("Model location: buffalo_pheno_mixed_corr_mvmt_10rs_shrubgrassref_wet_", Sys.Date(), ".rds")

Plot both models

# plot_model(buffalo.pheno.dry.mvmt)
# plot_model(buffalo.pheno.wet.mvmt)
# 
# plot_model(buffalo.pheno.dry.mvmt, transform = NULL)
# plot_model(buffalo.pheno.wet.mvmt, transform = NULL)

Create manuscript figures

Dry season

coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
                      "Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass", 
                                      "sl", "log(sl)", "cos(ta)", 
                                      "sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
                                      "log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
                                      "cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
                                      ),
                      "Estimate" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"],
                      "SE" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Std. Error"]
) %>%
  
  mutate(
    # Calculate confidence intervals
    LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE, 
    UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
    LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
    UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
    LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
    UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
    # Compute p-values
    pvalue = 2 * pnorm(-abs(Estimate) / SE)
  )

# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
  if (coef_df$pvalue[x] <= 0.001){
    return("***")
  } else if (coef_df$pvalue[x] <= 0.01) {
    return("**")
  } else if (coef_df$pvalue[x] <= 0.05) {
    return("*")
  }
})

# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")

# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))

# Specify the order in which the coefficients should be plotted
order <- c(
  "Floodplain",
  "DryGrassland",
  "OpenForest",
  "TreeGrass",
  "sl",
  "log(sl)",
  "cos(ta)",
  "sl:Floodplain",
  "sl:DryGrassland",
  "sl:OpenForest",
  "sl:TreeGrass",
  "log(sl):Floodplain",
  "log(sl):DryGrassland",
  "log(sl):OpenForest",
  "log(sl):TreeGrass",
  "cos(ta):Floodplain",
  "cos(ta):DryGrassland",
  "cos(ta):OpenForest",
  "cos(ta):TreeGrass"
)

# Specify the order in which the coefficients should be plotted
order_no_mvmt <- c(
  "Floodplain",
  "DryGrassland",
  "OpenForest",
  "TreeGrass"
)

# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

Plot the coefficients

min_x <- -0.5
max_x <- 0.8

# coefficients on the x-axis
dry_mvmt_plot <- ggplot(data = coef_df, 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 12.5, yend   = 12.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 15.5, yend   = 15.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  # drop the filter argument when plotting all covariates
  geom_errorbarh(data = coef_df2,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  scale_y_discrete(limits = rev(order)) + 
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#FF8123", "#9B4200")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  # ggtitle("Dry season") +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#FF8123")))

dry_mvmt_plot

Warning: Removed 9 rows containing missing values or values outside the scale range
(`geom_text()`).

ggsave(paste0("outputs/plots/pheno_model_mixed_corr_dry_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 150, units = "mm", scale = 1, dpi = 600)

Warning: Removed 9 rows containing missing values or values outside the scale range
(`geom_text()`).

Drop the filter argument when plotting all covariates

dry_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")), 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  # drop the filter argument when plotting all covariates
  geom_errorbarh(data = coef_df2_no_mvmt,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  # scale_y_discrete(limits = rev(order)) + 
  scale_y_discrete(limits = rev(order_no_mvmt)) + 
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#FF8123", "#9B4200")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#FF8123")))

dry_mvmt_plot

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

ggsave(paste0("outputs/plots/pheno_model_mixed_corr_dry_mvmt_WITHOUT_MVMT_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 60, units = "mm", scale = 1, dpi = 600)

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

Wet season

coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
                      "Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass", 
                                      "sl", "log(sl)", "cos(ta)", 
                                      "sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
                                      "log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
                                      "cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
                                      ),
                      "Estimate" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"],
                      "SE" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Std. Error"]
) %>%
  
  mutate(
    # Calculate confidence intervals
    LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE, 
    UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
    LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
    UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
    LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
    UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
    # Compute p-values
    pvalue = 2 * pnorm(-abs(Estimate) / SE)
  )

# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
  if (coef_df$pvalue[x] <= 0.001){
    return("***")
  } else if (coef_df$pvalue[x] <= 0.01) {
    return("**")
  } else if (coef_df$pvalue[x] <= 0.05) {
    return("*")
  }
})

# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")

# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))

# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

Plot the coefficients

# coefficients on the x-axis
wet_mvmt_plot <- ggplot(data = coef_df, 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 12.5, yend   = 12.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 15.5, yend   = 15.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  geom_point(shape = 3, size = 2.5) +
  geom_errorbarh(data = coef_df2,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  scale_y_discrete(limits = rev(order)) +
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#238DFF", "#004691")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#1377E2")))

wet_mvmt_plot

Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).

ggsave(paste0("outputs/plots/pheno_model_mixed_corr_wet_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 150, units = "mm", scale = 1, dpi = 600)

Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).

Drop the filter argument when plotting all covariates

wet_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")), 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  geom_errorbarh(data = coef_df2_no_mvmt,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  scale_y_discrete(limits = rev(order_no_mvmt)) + 
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#238DFF", "#004691")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#1377E2")))

wet_mvmt_plot

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

ggsave(paste0("outputs/plots/pheno_model_mixed_corr_wet_mvmt_WITHOUT_MVMT_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 60, units = "mm", scale = 1, dpi = 600)

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

Updating movement parameters, with vegetation interactions

# from the step generation script
tentative_shape <- 0.438167
tentative_scale <- 534.3507
tentative_kappa <- 0.1848126

sl_extent <- 1000 # metres
sl_dry <- data.frame(x = seq(0, sl_extent, length.out = 1001))

Dry season

coef_df_dry_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
                           "scaled_coefficient" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"])

# step length coefficients
sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate)) 
sl_dry_mvmt_mvparams <- sl_dry_mvmt_mvparams %>% mutate(
  scale_factor = sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
  
# log step length coefficients
log_sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate)) 
log_sl_dry_mvmt_mvparams <- log_sl_dry_mvmt_mvparams %>% mutate(
  scale_factor = log_sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# turning angle coefficients
cos_ta_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate)) 
cos_ta_dry_mvmt_mvparams <- cos_ta_dry_mvmt_mvparams %>% mutate(
  scale_factor = cos_ta_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

tentative_mean_sl <- tentative_shape * tentative_scale

# shrub_grass
shrub_grass_mvmt_dry <- data.frame(
  "pheno" = "ShrubGrass",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1],
  "updated_scale" = 1/((1/tentative_scale) - sl_dry_mvmt_mvparams$nat_coef[1]),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1]
)

shrub_grass_mean_sl_dry <- shrub_grass_mvmt_dry$updated_shape * shrub_grass_mvmt_dry$updated_scale

# floodplain
floodplain_mvmt_dry <- data.frame(
  "pheno" = "Floodplain",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[2],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[2])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[2]
)

floodplain_mean_sl_dry <- floodplain_mvmt_dry$updated_shape * floodplain_mvmt_dry$updated_scale

# dry_grassland
dry_grassland_mvmt_dry <- data.frame(
  "pheno" = "DryGrassland",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[3],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[3])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[3]
)

dry_grassland_mean_sl_dry <- dry_grassland_mvmt_dry$updated_shape * dry_grassland_mvmt_dry$updated_scale

# open_forest
open_forest_mvmt_dry <- data.frame(
  "pheno" = "OpenForest",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[4],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[4])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[4]
)

open_forest_mean_sl_dry <- open_forest_mvmt_dry$updated_shape * open_forest_mvmt_dry$updated_scale

# tree_grass
tree_grass_mvmt_dry <- data.frame(
  "pheno" = "TreeGrass",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[5],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[5])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[5]
)

tree_grass_mean_sl_dry <- tree_grass_mvmt_dry$updated_shape * tree_grass_mvmt_dry$updated_scale

Creating Gamma distributions

sl_dry_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))

# Tentative distribution
sl_dry_mvmt$Tentative <- dgamma(x = sl_dry_mvmt$x,
                                 shape = tentative_shape,
                                 scale = tentative_scale)
# Shrub grass
sl_dry_mvmt$ShrubGrass <- dgamma(x = sl_dry_mvmt$x,
                             shape = shrub_grass_mvmt_dry$updated_shape,
                             scale = shrub_grass_mvmt_dry$updated_scale)
# Bare ground
sl_dry_mvmt$Floodplain <- dgamma(x = sl_dry_mvmt$x,
                             shape = floodplain_mvmt_dry$updated_shape,
                             scale = floodplain_mvmt_dry$updated_scale)
# Dry grassland
sl_dry_mvmt$DryGrassland <- dgamma(x = sl_dry_mvmt$x,
                             shape = dry_grassland_mvmt_dry$updated_shape,
                             scale = dry_grassland_mvmt_dry$updated_scale)
# Open forest
sl_dry_mvmt$OpenForest <- dgamma(x = sl_dry_mvmt$x,
                             shape = open_forest_mvmt_dry$updated_shape,
                             scale = open_forest_mvmt_dry$updated_scale)
# Tree grass
sl_dry_mvmt$TreeGrass <- dgamma(x = sl_dry_mvmt$x,
                             shape = tree_grass_mvmt_dry$updated_shape,
                             scale = tree_grass_mvmt_dry$updated_scale)

# Pivot from wide to long data
sl_dry_mvmt <- sl_dry_mvmt %>% 
  pivot_longer(cols = -x)

# Specify the order in which the coefficients should be plotted
order <- c(
  "ShrubGrass",
  "Floodplain",
  "DryGrassland",
  "OpenForest",
  "TreeGrass",
  "Tentative"
)

sl_dry_mvmt$name <- factor(sl_dry_mvmt$name, levels = order)


mean_values <- c("ShrubGrass" = round(shrub_grass_mean_sl_dry, 1),
                 "Floodplain" = round(floodplain_mean_sl_dry, 1), 
                 "DryGrassland" = round(dry_grassland_mean_sl_dry, 1), 
                 "OpenForest" = round(open_forest_mean_sl_dry, 1), 
                 "TreeGrass" = round(tree_grass_mean_sl_dry, 1), 
                 "Tentative" = round(tentative_mean_sl, 1))

# Construct dynamic labels
labels <- paste(names(mean_values), ": Mean SL = ", mean_values, sep = "")

# Plot
ggplot(sl_dry_mvmt, aes(x = x, y = value, color = name, linetype = name)) +
  geom_line(alpha = 0.75, linewidth = 1) +
  scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
  scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black"),
                      labels = labels) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted"),
                        labels = labels) +
  theme_bw() +
  theme(legend.position = c(0.65, 0.65))

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

ggsave(paste0("outputs/plots/pheno_model_mixed_corr_sl_dry_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 1000)

Creating von Mises distributions

ta_dry_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))

# Tentative distribution
ta_dry_mvmt$Tentative <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = tentative_kappa)
# Shrub grass
ta_dry_mvmt$ShrubGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = shrub_grass_mvmt_dry$updated_kappa)
# Bare ground
ta_dry_mvmt$Floodplain <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = floodplain_mvmt_dry$updated_kappa)
# Dry grassland
ta_dry_mvmt$DryGrassland <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = dry_grassland_mvmt_dry$updated_kappa)
# Open forest
ta_dry_mvmt$OpenForest <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = open_forest_mvmt_dry$updated_kappa)
# Tree grass
ta_dry_mvmt$TreeGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = tree_grass_mvmt_dry$updated_kappa)

# Pivot from wide to long data
ta_dry_mvmt <- ta_dry_mvmt %>% 
  pivot_longer(cols = -x)

ta_dry_mvmt$name <- factor(ta_dry_mvmt$name, levels = order)

# Plot
ggplot(ta_dry_mvmt, 
       aes(x = x, y = value, color = factor(name), linetype = factor(name))
       ) +
  geom_line(alpha = 0.75, linewidth = 1) +
  xlab("Turning angle (radians)") +
  scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
  scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
                     labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black")) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted")) +
  # scale_colour_viridis_d("Habitat", option = "D") +
  # ggtitle("Dry season") +
  theme_bw() +
  theme(legend.position = "bottom")

ggsave(paste0("outputs/plots/pheno_model_mixed_corr_ta_dry_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 600)

Wet season

coef_df_wet_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
                                        "scaled_coefficient" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"])

# step length coefficients
sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate)) 
sl_wet_mvmt_mvparams <- sl_wet_mvmt_mvparams %>% mutate(
  scale_factor = sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# log step length coefficients
log_sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate)) 
log_sl_wet_mvmt_mvparams <- log_sl_wet_mvmt_mvparams %>% mutate(
  scale_factor = log_sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# turning angle coefficients
cos_ta_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate)) 
cos_ta_wet_mvmt_mvparams <- cos_ta_wet_mvmt_mvparams %>% mutate(
  scale_factor = cos_ta_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# shrub_grass
shrub_grass_mvmt_wet <- data.frame(
  "pheno" = "ShrubGrass",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1],
  "updated_scale" = 1/((1/tentative_scale) - sl_wet_mvmt_mvparams$nat_coef[1]),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1]
)

shrub_grass_mean_sl_wet <- shrub_grass_mvmt_wet$updated_shape * shrub_grass_mvmt_wet$updated_scale

# floodplain
floodplain_mvmt_wet <- data.frame(
  "pheno" = "Floodplain",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[2],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[2])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[2]
)

floodplain_mean_sl_wet <- floodplain_mvmt_wet$updated_shape * floodplain_mvmt_wet$updated_scale

# dry_grassland
dry_grassland_mvmt_wet <- data.frame(
  "pheno" = "wetGrassland",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[3],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[3])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[3]
)

dry_grassland_mean_sl_wet <- dry_grassland_mvmt_wet$updated_shape * dry_grassland_mvmt_wet$updated_scale

# open_forest
open_forest_mvmt_wet <- data.frame(
  "pheno" = "OpenForest",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[4],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[4])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[4]
)

open_forest_mean_sl_wet <- open_forest_mvmt_wet$updated_shape * open_forest_mvmt_wet$updated_scale

# tree_grass
tree_grass_mvmt_wet <- data.frame(
  "pheno" = "TreeGrass",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[5],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[5])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[5]
)

tree_grass_mean_sl_wet <- tree_grass_mvmt_wet$updated_shape * tree_grass_mvmt_wet$updated_scale

Creating Gamma distributions

sl_wet_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))

# Tentative distribution
sl_wet_mvmt$Tentative <- dgamma(x = sl_wet_mvmt$x,
                                shape = tentative_shape,
                                scale = tentative_scale)
# Shrub grass
sl_wet_mvmt$ShrubGrass <- dgamma(x = sl_wet_mvmt$x,
                                 shape = shrub_grass_mvmt_wet$updated_shape,
                                 scale = shrub_grass_mvmt_wet$updated_scale)
# Bare ground
sl_wet_mvmt$Floodplain <- dgamma(x = sl_wet_mvmt$x,
                                 shape = floodplain_mvmt_wet$updated_shape,
                                 scale = floodplain_mvmt_wet$updated_scale)
# Dry grassland
sl_wet_mvmt$DryGrassland <- dgamma(x = sl_wet_mvmt$x,
                                   shape = dry_grassland_mvmt_wet$updated_shape,
                                   scale = dry_grassland_mvmt_wet$updated_scale)
# Open forest
sl_wet_mvmt$OpenForest <- dgamma(x = sl_wet_mvmt$x,
                                 shape = open_forest_mvmt_wet$updated_shape,
                                 scale = open_forest_mvmt_wet$updated_scale)
# Tree grass
sl_wet_mvmt$TreeGrass <- dgamma(x = sl_wet_mvmt$x,
                                shape = tree_grass_mvmt_wet$updated_shape,
                                scale = tree_grass_mvmt_wet$updated_scale)

# Pivot from wide to long data
sl_wet_mvmt <- sl_wet_mvmt %>% 
  pivot_longer(cols = -x)

sl_wet_mvmt$name <- factor(sl_wet_mvmt$name, levels = order)

mean_values_wet <- c("ShrubGrass" = round(shrub_grass_mean_sl_wet, 1),
                     "Floodplain" = round(floodplain_mean_sl_wet, 1), 
                     "DryGrassland" = round(dry_grassland_mean_sl_wet, 1), 
                     "OpenForest" = round(open_forest_mean_sl_wet, 1), 
                     "TreeGrass" = round(tree_grass_mean_sl_wet, 1), 
                     "Tentative" = round(tentative_mean_sl, 1))

# Construct dynamic labels
labels <- paste(names(mean_values_wet), ": Mean SL = ", mean_values_wet, sep = "")

# Plot
ggplot(sl_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
  geom_line(alpha = 0.75, linewidth = 1) +
  scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
  scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black"),
                      labels = labels) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted"),
                        labels = labels) +
  # scale_colour_viridis_d("Habitat", option = "D") +
  # ggtitle("Wet season") +
  theme_bw() +
  theme(legend.position = c(0.65, 0.65))

ggsave(paste0("outputs/plots/pheno_model_mixed_corr_sl_wet_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 1000)

Creating von Mises distributions

ta_wet_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))

# Tentative distribution
ta_wet_mvmt$Tentative <- circular::dvonmises(x = ta_wet_mvmt$x,
                                             mu = 0,
                                             kappa = tentative_kappa)
# Shrub grass
ta_wet_mvmt$ShrubGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
                                              mu = 0,
                                              kappa = shrub_grass_mvmt_wet$updated_kappa)
# Bare ground
ta_wet_mvmt$Floodplain <- circular::dvonmises(x = ta_wet_mvmt$x,
                                              mu = 0,
                                              kappa = floodplain_mvmt_wet$updated_kappa)
# Dry grassland
ta_wet_mvmt$DryGrassland <- circular::dvonmises(x = ta_wet_mvmt$x,
                                                mu = 0,
                                                kappa = dry_grassland_mvmt_wet$updated_kappa)
# Open forest
ta_wet_mvmt$OpenForest <- circular::dvonmises(x = ta_wet_mvmt$x,
                                              mu = 0,
                                              kappa = open_forest_mvmt_wet$updated_kappa)
# Tree grass
ta_wet_mvmt$TreeGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
                                             mu = 0,
                                             kappa = tree_grass_mvmt_wet$updated_kappa)

# Pivot from wide to long data
ta_wet_mvmt <- ta_wet_mvmt %>% 
  pivot_longer(cols = -x)

ta_wet_mvmt$name <- factor(ta_wet_mvmt$name, levels = order)

# Plot
ggplot(ta_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
  geom_line(alpha = 0.75, linewidth = 1) +
  xlab("Turning angle (radians)") +
  scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
  scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
                     labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black")) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted")) +
  # scale_colour_viridis_d("Habitat", option = "D") +
  # ggtitle("Dry season") +
  theme_bw() +
  theme(legend.position = "bottom")

ggsave(paste0("outputs/plots/pheno_model_mixed_corr_ta_wet_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 600)

Session Information

sessionInfo()

R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Hobart
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] amt_0.2.2.0         circular_0.5-1      lemon_0.5.0        
 [4] beepr_2.0           tictoc_1.2.1        glmmTMB_1.1.11     
 [7] sjPlot_2.8.17       TwoStepCLogit_1.2.6 lubridate_1.9.4    
[10] forcats_1.0.0       stringr_1.5.1       dplyr_1.1.4        
[13] purrr_1.0.4         readr_2.1.5         tidyr_1.3.1        
[16] tibble_3.2.1        ggplot2_3.5.2       tidyverse_2.0.0    

loaded via a namespace (and not attached):
 [1] sjlabelled_1.2.0    tidyselect_1.2.1    farver_2.1.2       
 [4] fastmap_1.2.0       bayestestR_0.15.3   sjstats_0.19.0     
 [7] digest_0.6.37       timechange_0.3.0    lifecycle_1.0.4    
[10] sf_1.0-20           survival_3.8-3      magrittr_2.0.3     
[13] compiler_4.5.0      rlang_1.1.6         tools_4.5.0        
[16] yaml_2.3.10         knitr_1.50          labeling_0.4.3     
[19] bit_4.6.0           classInt_0.4-11     plyr_1.8.9         
[22] RColorBrewer_1.1-3  KernSmooth_2.23-26  withr_3.0.2        
[25] numDeriv_2016.8-1.1 grid_4.5.0          datawizard_1.1.0   
[28] e1071_1.7-16        scales_1.4.0        MASS_7.3-65        
[31] insight_1.2.0       cli_3.6.5           mvtnorm_1.3-3      
[34] crayon_1.5.3        rmarkdown_2.29      ragg_1.4.0         
[37] reformulas_0.4.1    generics_0.1.3      rstudioapi_0.17.1  
[40] performance_0.13.0  tzdb_0.5.0          parameters_0.25.0  
[43] minqa_1.2.8         DBI_1.2.3           proxy_0.4-27       
[46] audio_0.1-11        splines_4.5.0       parallel_4.5.0     
[49] effectsize_1.0.0    vctrs_0.6.5         boot_1.3-31        
[52] Matrix_1.7-3        jsonlite_2.0.0      hms_1.1.3          
[55] bit64_4.6.0-1       systemfonts_1.2.3   units_0.8-7        
[58] glue_1.8.0          nloptr_2.2.1        stringi_1.8.7      
[61] gtable_0.3.6        ggeffects_2.2.1     lme4_1.1-37        
[64] pillar_1.10.2       htmltools_0.5.8.1   R6_2.6.1           
[67] TMB_1.9.17          textshaping_1.0.1   Rdpack_2.6.4       
[70] vroom_1.6.5         evaluate_1.0.3      lattice_0.22-6     
[73] rbibutils_2.3       class_7.3-23        Rcpp_1.0.14        
[76] gridExtra_2.3       nlme_3.1-168        mgcv_1.9-1         
[79] xfun_0.52           sjmisc_2.8.10       pkgconfig_2.0.3